CRS
Diego Hernangómez
Un sistema de referencia de coordenadas (o CRS por sus siglas en inglés, Coordinate Reference System) permite relacionar datos espaciales con su localización en la superficie terrestre. Los CRS constituyen por tanto un aspecto fundamental en el análisis y representación de datos espaciales, ya que nos permiten identificar con exactitud la posición de los datos sobre el globo terráqueo.
Así mismo, cuando se trabaja con datos espaciales provenientes de distintas fuentes de información, es necesario comprobar que dichos datos se encuentran definidos en el mismo CRS:
En el ejemplo anterior, ambos puntos (verde y rojo) presentan los mismos valores de coordenadas en los ejes X e Y, en este caso las correspondientes a la ciudad de Toledo.
Sin embargo, presentan distintos CRS. Por este motivo, al representar ambos puntos en un mapa, se observa que no se están refiriendo a la misma localización geográfica. Esto es así porque el CRS define la referencia (punto x=0 e y =0) y las unidades de los ejes (grados, metros, millas).
Como conclusión, además de disponer de las coordenadas de los datos espaciales, es necesario conocer el CRS en el que están definidos para conocer de manera exacta su localización geográfica. Además, nótese que para cualquier análisis de datos espaciales es necesario que todos los geodatos se encuentren referenciados en el mismo CRS. Esto se consigue transformando (o proyectando) los datos a un CRS común, nunca sobreescribiendo el CRS de los mismos.
Tipos de CRS
A continuación se definen los dos grandes tipos de CRS, los CRS geográficos y los CRS proyectados.
CRS geográficos
Los CRS geográficos son aquellos en los que los parámetros empleados para localizar una posición espacial son la latitud y la longitud:
Latitud: Es la distancia angular expresada en grados sobre el plano definido por el ecuador terrestre. Determina la posición sobre de una localización en el eje Norte-Sur de la Tierra y toma valores en el rango \([-90º,90º]\) . Las líneas imaginarias determinadas por una sucesión de puntos con la misma latitud a lo largo del eje Este-Oeste se denominan paralelos:
Longitud: Es la distancia angular expresada en grados sobre el plano definido por el meridiano de Greenwich. Determina la posición sobre de una localización en el eje Este-Oeste de la Tierra y toma valores en el rango \([-180º,180º]\) . Las líneas imaginarias determinadas por una sucesión de puntos con la misma longitud a lo largo del eje Este-Oeste se denominan meridianos:
Es muy importante destacar que en un sistema de coordenadas geográfico, es decir, basado en latitudes y longitudes, las distancias entre dos puntos representan distancias angulares. Por ejemplo, la distancia entre el meridiano de Greenwich y el meridiano correspondiente a la longitud 20º siempre es de +20º. Sin embargo, debido a la forma esférica de la Tierra, la longitud en metros entre ambos meridianos no es constante:
CRS proyectados
La representación de formas tridimensionales en un soporte plano (dos dimensiones) presenta algunos retos. Por ello, es habitual trabajar con proyecciones de mapas.
Una proyección geográfica es un método para reducir la superficie de la esfera terrestre a un sistema cartesiano de dos dimensiones. Para ello, es necesario transformar las coordenadas longitud y latitud en coordenadas cartesianas x e y.
Es importante destacar que las proyecciones pueden incluir un punto de origen (X=0, Y=0) y unas unidades de distancia (habitualmente metros) específicas. Por ejemplo, la proyección cónica equiáreas de Albers (específica para Estados Unidos) define su punto de referencia (0,0) en la latitud 40º N y longitud 96º, y la unidad de variación están definida en metros. De ahí la importancia de conocer el CRS de los datos geográficos, como se expuso al principio de este tema.
Existen varias familias de proyecciones, que se pueden clasificar de diversas maneras:
Por tipo de superficie de proyección
El proceso de trasladar puntos de una esfera a un plano puede plantearse de manera práctica como el ejercicio de envolver una esfera con una superificie plana (como una hoja de papel) y trasladar los puntos de la esfera de manera lineal al punto de la superficie plana más cercano a ella.
A partir de este ejercicio, se plantean tres posibles soluciones, dependiendo del tipo de superficie que se use para proyectar:
- Proyecciones cilíndricas: Son aquellas proyecciones donde la superficie de proyección conforma un cilindro alrededor de la Tierra. Una de las proyecciones cilíndricas más conocidas es la proyección de Mercator.
Proyecciones cónicas: En este tipo de proyecciones, se plantea la superficie de proyección como una forma cónica. Como ejemplo, la proyección cónica equiáreas de Albers es una de las proyecciones que más suele usarse en la representación de mapas de América del Norte:
Proyecciones acimutales o planares: En este tipo de proyección se proyecta una porción de la Tierra sobre un plano que es tangente a la misma en el punto de referencia. Como ejemplos de proyecciones acimutales podemos destacar la proyección ortográfica:
Por métrica a preservar
Es importante tener en cuenta que cualquier proyección de la superficie de la Tierra produce distorsiones en una o varias características geográficas. Como ejemplos clásicos, la proyección de Mercator produce distorsiones del área especialmente en aquellas regiones más cercanas a los polos (Groenlandia, que la proyección de Mercator presenta una área similar a la de África, presenta menor superificie real que Argelia). Otras de las métricas que suele verse distorsionada son la distancia entre dos puntos geográficos, la dirección o la forma de regiones de la Tierra.
A lo largo de la Historia se han desarrollado diversas proyecciones cuyo objetivo es preservar alguna o varias de las propiedades mencionadas anteriormente, sin embargo es importante destacar que no existe una proyección que sea capaz de preservar todas las métricas a la vez.
Según la metrica a presevar, las proyecciones se pueden clasificar en:
Proyecciones conformales: Estas proyecciones intentan preservar los ángulos que se forman en la superficie terrestre. Por ejemplo, la proyección de Mercator representa ángulos rectos en las intersecciones de los paralelos y los meridianos.
Proyecciones equivalentes: Estas proyecciones preservan las proporciones de las áreas, provocando a su vez deformaciones en el resto de características, como la forma o los ángulos. La proyección acimutal equivalente de Lambers es un tipo de proyección equivalente.
Proyecciones equidistantes: Este tipo de proyección preserva la distancia entre dos puntos geográficos específicos. Por ejemplo, la proyección Plate carré preserva la distancia entre el Polo Norte y el Polo Sur.
- Proyecciones de compromiso: Este tipo de proyección no intenta preservar ninguna métrica en concreto. En su lugar, se centran en intentar encontrar un equilibrio entre las diversas distorsiones que provocan para intentar dar una representación más o menos representativa de la superficie terrestre. La proyección de Winkel Tripel, usada en los mapas de National Geographic, es un ejemplo de proyección de compromiso.
En los anteriores ejemplos se ha añadido a cada proyección la indicatriz de Tissot. Consiste en una serie de círculos imaginarios de igual área distribuidos sobre la superficie esférica de la Tierra en determinados puntos. De este manera, al presentar la indicatriz de Tissot en una proyección específica, se puede entender de una manera intuitiva la distorsión provocada por dicha proyección, ya que los círculos se ven distorsionados o preservados según los parámetros y la naturaleza de la proyección en cuestión.
Trabajando con proyecciones en R
Existe toda una serie de proyecciones predefinidas, identificadas mediante los códigos EPSG, ESRI, WKT o proj4 (en desuso en R, pero todavía admitidos). Existen varios recursos web donde se pueden consultar y seleccionar los códigos correspondientes:
El paquete sf permite obtener los parámetros de estas proyecciones mediante la función st_crs():
library(sf)
# Ejemplo: EPSG WGS 84 (Sistema Global GPS): EPSG 4326
st_crs(4326)## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["unknown"],
## AREA["World"],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]# Usando código ESRI North America Albers Equal Area Conic
st_crs("ESRI:102008")## Coordinate Reference System:
## User input: ESRI:102008
## wkt:
## PROJCRS["North_America_Albers_Equal_Area_Conic",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["North_America_Albers_Equal_Area_Conic",
## METHOD["Albers Equal Area",
## ID["EPSG",9822]],
## PARAMETER["Latitude of false origin",40,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-96,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",20,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",60,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["North America - Canada and USA (CONUS, Alaska mainland)"],
## BBOX[23.81,-172.54,86.46,-47.74]],
## ID["ESRI",102008]]# Usando proj4string: Robinson
st_crs("+proj=robin")## Coordinate Reference System:
## User input: +proj=robin
## wkt:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["unknown",
## METHOD["Robinson"],
## PARAMETER["Longitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]De esta manera, es posible proyectar un objeto sf mediante la función st_transform():
# Usa datos del paquete mapSpain
library(giscoR)
paises <- gisco_get_countries()
# Comprobamos el CRS de estos datos
# Se puede almacenar en un objeto y usar posteriormente
st_crs(paises)## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCS["WGS 84",
## DATUM["WGS_1984",
## SPHEROID["WGS 84",6378137,298.257223563,
## AUTHORITY["EPSG","7030"]],
## AUTHORITY["EPSG","6326"]],
## PRIMEM["Greenwich",0,
## AUTHORITY["EPSG","8901"]],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4326"]]# Plot
plot(st_geometry(paises), axes = TRUE)
# Proyectamos a Mercator
# El eje cambia porque Mercator usa metros
paises_merc <- st_transform(paises, st_crs(3857))
plot(st_geometry(paises_merc), axes = TRUE)
# Proyectamos a Robinson
paises_robin <- st_transform(paises, st_crs("+proj=robin"))
plot(st_geometry(paises_robin), axes = TRUE)
Como se comentó anteriormente, cuando se usan geodatos de diversas fuentes, es necesario que todos presenten el mismo CRS. En este ejemplo se muestra lo que ocurre si esto no se cumple:
# Añadimos a este mapa puertos mundiales de giscoR
puertos <- gisco_get_ports()
plot(st_geometry(paises_robin), main = "Puertos en el mundo")
plot(st_geometry(puertos), add = TRUE, col = "red", pch = 20)
# Ha habido algun error... Comprueba CRS
st_crs(puertos) == st_crs(paises_robin)## [1] FALSE# Los puertos no están en Robinson! Proyectamos al mismo CRS
puertos_robin <- st_transform(puertos, st_crs(paises_robin))
plot(st_geometry(paises_robin), main = "Puertos en el mundo")
plot(st_geometry(puertos_robin), add = TRUE, col = "blue", pch = 20)
Como vemos, en el primer mapa los puertos se concentran en un único punto, dado que no están referenciados en el mismo CRS. Tras proyectarlos, el mapa se representa adecuadamente.
En otros paquetes, como sp o raster, existen funciones parecidas. Cuando empleemos el paquete sp podemos usar las funciones CRS() y spTransform():
library(sp)
# Convertimos sf a sp
paises_sp <- as(paises, "Spatial")
paises_sp@proj4string## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCS["WGS 84",
## DATUM["WGS_1984",
## SPHEROID["WGS 84",6378137,298.257223563,
## AUTHORITY["EPSG","7030"]],
## AUTHORITY["EPSG","6326"]],
## PRIMEM["Greenwich",0,
## AUTHORITY["EPSG","8901"]],
## UNIT["degree",0.0174532925199433,
## AUTHORITY["EPSG","9122"]],
## AUTHORITY["EPSG","4326"]]# En sp podemos usar:
CRS("+proj=robin")## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["unknown",
## METHOD["Robinson"],
## PARAMETER["Longitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]# O también desde sf
CRS(st_crs(paises_robin)$proj4string)## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["unknown",
## METHOD["Robinson"],
## PARAMETER["Longitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]paises_sp_robin <- spTransform(paises_sp, CRS("+proj=robin"))
plot(paises_sp_robin)
En el caso de un objeto raster, podemos usar crs() y projectRaster():
library(raster)
# Extrae información de altitud para Paises Bajos
elev <- getData("alt", country = "NLD", path = tempdir())
raster::crs(elev)## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]plot(elev)
# Transforma
elev_robinson <- projectRaster(elev, crs = crs("+proj=robin"))
plot(elev_robinson)
Por último, en el paquete terra las funciones correspondientes son crs() y project():
library(terra)
# Convierte de raster a terra
elev_terra <- rast(elev)
cat(terra::crs(elev_terra))## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]plot(elev_terra)
# Transforma
elev_terra_robinson <- terra::project(elev_terra, terra::crs(elev_terra))
plot(elev_terra_robinson)
¿Qué proyección uso?
El CRS adecuado para cada análisis depende de la localización y el rango espacial de los datos. Un CRS adecuado para representar un mapa del mundo puede no serlo para representar datos de zonas específicas de la Tierra. Los recursos web mencionados anteriormente permiten la búsqueda de CRS por zona geográfica, y adicionalmente en R existe el paquete crsuggest (Walker 2021) que nos facilita la labor, sugiriendo el CRS más adecuado para cada zona:
library(crsuggest)
# CRS para Paises Bajos
# Usando raster
suggest_crs(elev)## # A tibble: 10 × 6
## crs_code crs_name crs_type crs_gcs crs_units crs_proj4
## <chr> <chr> <chr> <dbl> <chr> <chr>
## 1 28992 Amersfoort / RD… project… 4289 m +proj=sterea +lat_0=52.…
## 2 28991 Amersfoort / RD… project… 4289 m +proj=sterea +lat_0=52.…
## 3 5651 ETRS89 / UTM zo… project… 4258 m +proj=tmerc +lat_0=0 +l…
## 4 5649 ETRS89 / UTM zo… project… 4258 m +proj=tmerc +lat_0=0 +l…
## 5 3812 ETRS89 / Belgia… project… 4258 m +proj=lcc +lat_0=50.797…
## 6 3447 ETRS89 / Belgia… project… 4258 m +proj=lcc +lat_0=50.797…
## 7 31370 Belge 1972 / Be… project… 4313 m +proj=lcc +lat_0=90 +lo…
## 8 31300 Belge 1972 / Be… project… 4313 m +proj=lcc +lat_0=90 +lo…
## 9 21500 Belge 1950 (Bru… project… 4809 m +proj=lcc +lat_0=90 +lo…
## 10 23095 ED50 / TM 5 NE project… 4230 m +proj=tmerc +lat_0=0 +l…# Probamos sugerencia
crs_suggest <- suggest_crs(elev, limit = 1)
elev_suggest <- projectRaster(elev, crs = raster::crs(crs_suggest$crs_proj4))
plot(elev_suggest)
# Ejemplo con sf: China
china <- gisco_get_countries(country = "China")
plot(st_geometry(china), axes = TRUE)
china_crs <- suggest_crs(china, limit = 1)
china_suggest <- st_transform(
china,
st_crs(as.integer(china_crs$crs_code))
)
plot(st_geometry(china_suggest), axes = TRUE)